R Notebook

Libraries

Load data

I am using the AHA dataset for writing my code because it is the simplest dataset I can work with. Simplest = significant number of baseline samples without issue of repeated measures, balanced distribution of classes, likely difference to be detected present.

p_ps = readRDS("/Users/aa370/Library/CloudStorage/Box-Box/project_davidlab/LAD_LAB_Personnel/Ammara_A/Projects/POMMS/R24_POMMS/data/processed/20221104_pomms_metadata_plant.rds")

p_ps

sample_data(p_ps)

Preprocess

Remove problematic samples

Taking out samples with less than or equal to 0 reads and keeping only entry timepoint reduces samples from 384 down to 240.

p_ps = p_ps %>%
  subset_samples( reads > 0) %>%
  subset_samples(timepoint == "Entry")

p_ps

Class distribution in this dataset: Case = 193 Control = 47

sample_data(p_ps) %>%
  as.data.frame() %>%
  as_tibble() %>%
  group_by(group) %>%
  summarise(count = n())

Data transform

otu_clr = abundances(p_ps, "clr") %>% # rclr transform could not be performed, get Nans
  t()

p_ps = phyloseq(otu_table(otu_clr, taxa_are_rows = FALSE),
                   sample_data(sample_data(p_ps)),
                   tax_table(tax_table(p_ps)))

rm(otu_clr)

Feature table

Extract data

Order of asv colnames in OTU table and tax table is the same.

features = p_ps@otu_table %>%
  as.data.frame() %>%
  as_tibble(rownames = NA) %>%
  rownames_to_column("sample")

feature_labels = c("sample")

 names = tax_table(p_ps) %>%
      data.frame() %>%
      lowest_level() %>%
   pull(name)

 feature_labels = append(feature_labels, names)

 feature_labels = make.unique(feature_labels) %>%
   make.names()

 colnames(features) = feature_labels

features

Attach data labels

#Subsampling: SMOTE

Zero variance features

No features have zero variance

nzv = subsampled %>%
  select(!group) %>%
  nearZeroVar(saveMetrics = TRUE)

rm(nzv)

Correlated predictors

summary(feature_cor[upper.tri(feature_cor)])

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.52293 -0.04691  0.10074  0.20495  0.35340  1.00000 

high_cor

  [1]  89  97 102 105 119 139 141 148 151 153 155 156 176 185 186 207 208 216 220 221 222 225 231 236
 [25] 238 241 251  32  55  69  75  76  77  78  80  83  85  86  88  91  92  94  95  98 101 103 104 106
 [49] 107 108 109 110 115 117 118 125 126 130 131 132 133 135 138 140 142 143 146 149 150 154 157 158
 [73] 159 100 160 172 174 179  74 181 189 190 191 192 193 194 195 199 200 202 203 204 205 206 209 210
 [97] 212 213 214 215 218 219 223 224 226 227 228 229 230 232 233 235 240 242 244 245 247 249

rm(high_cor)
rm(feature_cor)

Linear dependencies

feature_ld = filtered_subsampled %>%
  select(!group) %>%
  findLinearCombos()

filtered_subsampled = filtered_subsampled[, - feature_ld$remove]

Input table

Preprocessing reduces column number to 102 from 235

Modeling

Load functions from classification_function_AA file

results

[[1]]

[[2]]

NA

rm(label_list)

Plotting

AUC plot and t-test

auc_df %>%
  pivot_longer(cols = c("group", "group_shuffle"), names_to = "feature", values_to = "AUC") %>%
  ggplot()+
  geom_boxplot(aes(x=feature, y = AUC))+
  theme_classic()

t.test(auc_df$group, auc_df$group_shuffle)

    Welch Two Sample t-test

data:  auc_df$group and auc_df$group_shuffle
t = 21.57, df = 15.18, p-value = 8.37e-13
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.3677251 0.4482749
sample estimates:
mean of x mean of y 
    0.942     0.534 

Importance plots

look-up table

#Making asv reference for names in model
lookup_asv = tax_table(p_ps) %>% 
     data.frame() %>%
      lowest_level() %>%
   rownames_to_column("asv") %>% 
   select(asv, name) %>%
   mutate(name = make.names(make.unique(name)))

summarize importances

top 10 ranked list

# Getting the ranked list
ranking = importance_summary %>%
  filter(variable == "group_shuffle") %>%
  pull(food)

top10 = ranking[1:10]

Importance plot

importance_df %>%
  pivot_longer(cols = c("Helianthus.annuus":"Rheum.rhaponticum"), names_to = "food", values_to = "importance")  %>%
  filter(variable == "group_shuffle" & food %in% top10) %>%
  ggplot()+
  geom_boxplot(aes(x= reorder(food,-importance), y = importance))+
  labs(x= "Food", y ="Importance") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust=1))

Some of the important NAs for predictions only map to bacterial species in BLAST.

Important taxa directions

lookup_asv %>%
  filter(name == "Theobroma.cacao") %>%
  pull(asv)

[1] "ATCCTATTATTTTATTATTTTACGAAACTAAACAAAGGTTCAGCAAGCGAGAATAATAAAAAAAG"